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ABSTRACT 

The discovery of Jupiter-mass planets in close orbits about their parent stars has challenged models of planet 
formation. Recent observations have shown that a number of these planets have highly inclined, sometimes 
retrograde orbits about their parent stars, prompting much speculation as to their origin. It is known that 
migration alone cannot account for the observed population of these misaligned hot Jupiters, which suggests 
that dynamical processes after the gas disc dissipates play a substantial role in yielding the observed inclination 
and eccentricity distributions. One particularly promising candidate is planet-planet scattering, which is not 
very well understood in the non-linear regime of tides. Through three-dimensional hydrodynamical simulations 
of multi-orbit encounters, we show that planets that are scattered into an orbit about their parent stars with 
closest approach distance being less than approximately three times the tidal radius are either destroyed or 
completely ejected from the system. We find that as few as 5 and as many as 18 of the currently known 
hot Jupiters have a maximum initial apastron for scattering that lies well within the ice line, implying that 
these planets must have migrated either before or after the scattering event that brought them to their current 
positions. If stellar tides are unimportant (Q^ > 10^), disk migration is required to explain the existence of the 
hot Jupiters present in these systems. Additionally, we find that the disruption and/or ejection of Jupiter-mass 
planets deposits a Sun's worth of angular momentum onto the host star. For systems in which planet-planet 
scattering is common, we predict that planetary hosts have up to a 35% chance of possessing an obliquity 
relative to the invariable plane of greater than 90 degrees. 

Subject headings: Planet- star Interactions — Hydrodynamics — Gravitation — Chaos — Stars: rotation — 
Ultraviolet: planetary systems 



1. INTRODUCTION 

The search for planets about other stars has led to the dis- 
covery of dozens of planets with unusual properties. As both 
radial velocity and transit surveys are biased towards planets 
that are both massive and close to their parent stars, the re- 
gion of parameter space corresponding to planets with a mass 
larger than that of Neptune and a semi-major axis < 0.1 au 
is par t icularly well-explored ([Ida & Lin|2004[ |Shen & Turner 



2008 [ |Zakamska et al.j|2010| ). This has led to the discovery 
of many giant planets known colloquially as hot Jupiters and 
Neptunes, which are thought to have formed far away from 
their parent stars but were then later transplanted to their ob- 
served positions by currently undetermined means. Many of 
these exoplanets come so close to their parent stars that they 
toe the line between destruction and survival, with some ob- 
served exoplane ts in danger of being destroyed on a relatively 
short timescale ([Li et al.|2010| . Additionally, the inclination 
distribution of the hot Jupiters seems to demonstrate signif- 
icant mis alignment between the planet's orbit a nd the stellar 
spin axis ( [Triaud et al.|2010[ |Schlaufman|2010[ ), a surprising 
result that may require a dynamical process that acts after the 
protoplanetary disk dissipates. 

There are three primary physical processes that can deposit 
a planet on an orbit that is very close to its parent star: disk 
migration, the Kozai mechanism, and planet-planet scatter- 
ing. Disk migration can yield hot Jupiters, but as the star 
collapses from the same cloud as the protoplanetary disk that 
encircles it, it is difficult to explain the observed orbit mis- 
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alignments using this mechan ism alone (though see |Foucart & 
'Lai'2010VWats on et al.|2 010, for discussions regarding star- 
disk interactions). The Kozai mechanism (Kozai 1962) can 
lead to the large eccentricities required to produce close-in 
planets, but it can only operate in systems in which a massive 
planet or secondary star are present, and the mechanism may 
be mitigated by general relativistic effects that become im- 
portant before ti dal dissi pation is large enough to circularize 
the orbit (Taked a & Rasio 2005 ; Fabry cky & Tremaine 200^. 
Planet-planet scattering can produce both the observed semi- 
major axis and inclination distributions, and can deposit plan- 
ets close enough such that tides can circularize the orbits in a 
time that is less than the system age. Additionally, the object 
that acts as a scatterer can have approximately the same mass 
as the scattered object itself a nd still yield a hot Ju piter in a 
significant fraction of systems ( [Ford & Rasio|2008] ), negating 
the need for a non-planetary companion in the system. 

Previous hydrodynami cal work has only focused on the 
planet's first close fly-by ( |Faber et al.|2005[ hereafter FRW), 
and does not investigate how prolonged tidal forcing over 
many orbits affects a planet's chances for survival. In this pa- 
per we have performed hydrodynamical simulations of mul- 
tiple passages of a Jupiter-like planet by a Sun-like star, 
bridging the gap between numerical and analytical work 
that have focused on extremely close and extremely graz- 
ing encounters respectively. We find that scattering plan- 
ets into star-grazing orbits is more destructive than previ- 
ously thought, with Jupiter-like planets being destroyed or 
ejected at distances no smaller than 2.7 times the tidal ra- 
dius Tt = /?p(M*/Mp)^/^. As some exoplanets are currently 
observed to have semi-major axes less than twice this criti- 
cal value, their initial eccentricities may be required to have 
been substantially smaller than unity if planet-planet scatter- 
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ing is the mechanism responsible for bringing them so close 
to their host stars. This strongly suggests that planet-planet 
scattering alone cannot explain the complete observed popu- 
lation of close-in Jupiter-like exoplanets, and that the process 
must operate along with one of either the Kozai mechanism, 
disk migration, or both. These three processes likely act in 
concert to produce the observed population of hot gas giants, 
with the relative importance of each process being a function 
of the system's initial conditions. 

If planet-planet scattering is common enough to explain the 
existence of hot Jupiters, we predict that there should be two 
signatures of disruption that are readily detectable with to- 
day's instruments. Firstly, we find that the parent star can 
have its spin significantly altered by the accretion of material 
removed from the planet as a result of the disruption, produc- 
ing a star that can be significantly misaligned relative to any 
remaining planets. Secondly, we find that most planet disrup- 
tion events lead to the planet's ejection from the host system 
prior to the planet being completely destroyed, and that this 
ejected planet can remain almost as bright as its host star for 
centuries. 

In this paper we focus on the results of numerical hydrody- 
namical simulations that have been used to attempt to ascer- 
tain the true radii of destruction and ejection for Jupiter-like 
exoplanets, and the consequences of these planet-removing 
processes on their stellar hosts. In Section [2] we review the 
history of the analytical and numerical work done to char- 
acterize the orbital evolution of a planet that comes within 
a few tidal radii of its host star, and then we detail our par- 
ticular numerical approach to modeling tidal disruption. We 
report the results of our simulations in Section [3] In Section]?] 
we discuss the implications of our results, witn special atten- 
tion paid to the viability of various mechanisms for producing 
hot Jupiters, and the observational signatures of planetary dis- 
ruption and ejection. We summarize the shortcomings of our 
models and the possible fates of a Jupiter-like exoplanet in 
Section [5] Appendix |A| is provided to detail our algorithm 
used to simulate multiple orbits and for presenting tests of the 
algorithm's conservative properties. 

2. MODELING PLANETARY DISRUPTION 

2.1. Previous Disruption Models 

Tidal encounters between a point mass and an extended, 
initially spherical object have been described with progres- 
sively more detailed analytical models and simulations. The 
first analytica l models of tidal dissipation were laid out in 
Press & Teukolsky ( 1977 ) (hereafter PT), which assumes that 
the tidally excited body retains its spherical shape and that 
the motions induced by the encounter can be described with 
spherical harmonics. This model works quite well for encoun- 
ters where tides are weak, i.e. more than a few tidal radii away 
from the perturbing body, where the assumption of sphericity 
is valid. If mode-mode coupling is weak, the initially excited 
modes have no mechanism to share energy with each other 
during the encounter itself, and the PT formalism can accu- 
rately predict the amount of energy injected into the extended 
object. 

In the IPTI formalism, the amount of energy and angular 
momentum injected into an object is derived by taking the 
convolution of the object's modes of oscillation and the fre- 
quency decomposition of the tidal field, which drops the time- 
dependence from the governing equations. However, if the 
incoming object is already oscillating, the phase of these os- 



cillations relative to the time of periastron is important, which 
by construction cannot be resolved by the time-independent 
PT formalism. The absence of time-dependence also requires 
that the base- state is cylindrically symmetric relative to the 
line connecting the extended object to the perturber, which 
can only be accomplished for an initially non-rotating object 
by using a spherical geometry. As the amplitude of the mode 
oscillations approaches the size of the object itself, t he b ody 
must become highly non-spherical and therefore the |PT| for- 
malism must break down. 

These two shortcomings lead to the development of the 
"affine" model for tidal encounters ( [Carter & Luminet|1983| 
,1985) , which allows the initially spherical object to deform 
into a triaxial ellipsoid during and after the encounter. One 
axis of the ellipsoid is fixed to be perpendicular to the orbital 
plane, with the other two axes lying within the plane at a right 
angle to each other and with arbitrary orientation angle rela- 
tive to the first axis. This feature allows for the object to be 
followed dynamically, which means that the incoming state 
of the object's oscillations can affect the object's response to 
future encounters. 

Because the excitation of normal modes in a dynamically 
inert object can only yield an increase in the object's energy 
and angular momentum budgets, an initially spherical object 
would always find itself in a more-tightly bound orbit after the 
encounter. An important facet of the problem that the affine 
model can investigate is how the fundamental modes excited 
in previous encounters are de-excited in future encounters, 
which can lead to a positive change in the orbital energy. The 
inclusion of de-excitation of modes adds a chaotic compo- 
nent to the problem, with the orbital evolution behaving as a 
"random- walk" process ( Kochanek| 1 992^ |Mardling| 1 995] ) . 
While mass loss has been included ina nested-ellipsoid 



version of the affine model presented by Ivanov & Novikov 
(2001 ), asymmetrical mass loss is not treated as the mod 
els assume mirror symmetry. Asymmetric mass loss is ex- 
pected in real disruptions as the tidal field is stronger on the 
side of the object facing the point mass. This is due to the 
steep power-law dependence of the strength of the tidal field 
(oc r~^) and the non-linear evolution of the tide raised on the 
object, with both the velocity and amount of mass lost being 
larger on the side of th e objec t closest to the point mass. Nu- 
merical simulations by |FRW| show that this asymmetric mass 
loss leads to a substantial deviation from the nested-ellipsoid 
treatment, especially when the closest approach distance is 
< 2rt. Within this distance, the asymmetric removal of mass 
can lead to a positive increase in the orbital energy, which can 
cause the object to be completely ejected from the system if 

The disruption of polytropes in highly eccentric encoun- 
ters has bee n investigated numerically by both Lagrangian 
dNolthenius & Katz 1982; Bicknell & Gingold 1983; Evans] 
& Kochanek 1989; Kobayashi et al. 2004; Faber et al. 2005; 
Rossw og et al 2()08,.2009,.^^^ & Rosswog 2009i 

Lodato 'et al.||2009|) and Eulerian methods (Khokhlov e t"aL 
1993a b; 'Frolov et aLl[T994l |Diener et al.]^1997^ iGuillochon 
et al. 2009 ), with the principle focus being on stars or compact 
objects that are disrupted by point-like gravitational sources. 
However, most hydrodynamic simulations that have focused 
on the long-term survival of these systems only consider when 



the two objects have a mass ratio close to unity (Lee et al. 
2010[ |Loren-Aguilar et aT]|2009| ), or for non-disrupt i ve en 



counters ( Rathore 2005 1 ) . Recently, |Antonini et a 
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Figure 1. Final orbital energy £^orb scaled to the initial orbital energy £^orb,o (left panel), change in orbital energy A£^orb = ^orb,o~^orb scaled by £^orb,o (left 
sub-panel), and spin angular momentum /spin scaled by the characteristic angular momentum of Jupiter jj = GM^Rj (right panel) as functions of periastron 
distance rp after a single near-parabolic encounter between a Mp = Mj Jupiter-like planet and a M* = 10^ Mj star. The solid lines show the results from this work, 
with decreasing thickness corresponding to increasing maximum refinement, square markers show th e results of FRW and the dashed lines shows the prediction 
of PT using FRW s analytical n = 1 solutions for 7] and Si. Note that our results and that of |FRW| are consistent for rp < 2.5rt. At the highest resolution, our 
results are consistent with that of.PTiuntil rp > 4rt. 
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Figure 2. The double rainbow curves show orbital trajectories for both the planet (left) and star (right) during the encounter for different values of the periastron 
distance rp. The star's trajectory is magnified by a factor of M^/Mp = 10^ to make its motion apparent. While the orbits do precess slightly over the course of 
the simulations, this plot shows the orbits with the precession removed. Note that the r^ = 21 run (light green) experiences a particularly ejective encounter on its 
2nd periastron passage. 



performed low-resolution multiple-passage simulations in the 
context of the galactic center, allowing the exploration of a 
large parameter space at the expense of accuracy. 

To summarize, most analytical models of tides in an astro- 
physical context focus on objects which do not lose any mass 
at their closest approach, and conversely most numerical work 
has focused on encounters where the object is completely de- 
stroyed. The intermediate regime, where objects lose some 
mass but are not completely destroyed in their first passage, 
is largely uncharacterized. If the planet survives the initial 
encounter, some of the mass that is removed from the planet 
can become bound to the planet again as it recedes from pe- 
riastron to a region with a weaker tidal field. The return and 
subsequent re-accretion of this material is not treat ed at all 
in analytical models. As we will describe in Section 3.2 the 
structure of the planet can be significantly altered by the re- 
accretion of the planetary envelope, and the inclusion of this 
re-accretion into tidal disruption theory is necessary to deter- 
mine if a planet will ultimately survive. 

2.2. Our Approach 

While Lagrangian codes are well-suited for treating prob- 
lems where a hydrostatic object moves rapidly with respect to 



a fixed reference frame, their relatively poor scalability makes 
it difficult to follow the evolution of an object for many dy- 
namical timescales with sufficient spatial resolution. On the 
other hand, maintaining near-hydrostatic balance in rapidly 
advecting fra mes using Euleri an methods can also be quite 
challenging. [Robertson et a l.| (12009 ) show that performing 
Eulerian simulations in a boosted frame with unresolved pres- 
sure gradients leads to a spurious viscosity term that tends 
to damp out instabilities. This has been colloquially referred 
to as the non- Galilean invariance (GI) of the Rie mann prob- 



Tasker et 



aD [2QQ8l ). While [Robertson| 
iiTbe avoided with increased 



lem (SpringeT|[2QTQ| 
let al.l showed that GI issues can 
resolution, the requisite spatial resolution can be impractical 
for three-dimensional simulations that are evolved over thou- 
sands of dynamical timescales, or for when the bulk velocities 
are many times larger than the internal sound speeds within 
the simulation. The problem is compounded in the outermost 
layers of Jupiter-like planets, which exhibit particularly steep 
pressure gradients. 

To side-step the GI issues, our simulations are performed 
in the rest-frame of the planet, which limits the typical fluid 
velocities to the planet's sound speed for tidal disruption cal- 
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culations ^Guillochon et al.|r2QQ9 ). We model the planet as 
mi n= 1 , 1=2 polytrope, with the fluid being described by 
a polytropic equation of state (P ex p^), where the adiabatic 
index 7 = P. This gives a reasonable a pproximation to the 
interior structure of Jupiter-like planets ( Hubbard" 19 84 ), as 
long as the core is not a significant fraction of the planet's 
mass. Our simulations are constructed within the framework 



of FLASH ( [Fryxell et al.|2Q0 Q ), an adaptive-mesh, grid-based 
hydrodynamics code. We treat the star as a point-mass be- 
cause the distortion of the star itself contributes negligibly 
to the planet's orbital evolution in the case where Mp 
(iMatsumura et al. 2008). The relative positions of the star and 
tne planet are explicitly tracked over the course of the simula- 
tion using two virtual particles Xs (star) and Xp (planet), which 
are evolved alon gside the hy dro calculation using a Burlisch- 
Stoer integrator ( [Press et aL 1986) with maximum error con- 
strained to be less than 10"^-^. The planet's self-gravity is 
calculated using a multipole OQ^^^N) expansion of the fluid, 
where we take /max = 12. 

The flux-calculation step in any hydro code depends on 
the time-dependent force applied to the fluid over the dura- 
tion of the step. Because the gravitational field is applied 
as a source term in most multi-dimensional hydro codes that 
include self-gravity, conservation of energy is not explicitly 
achievable when self-gravity is included, as the potential at 
timestep m + 1 is unknown and must be estimated. In the case 
of self-bound objects in hydrostatic equilibrium, small per- 
turbations can lead to a systematic drift in the total energy 
of the system. For hydrostatic objects, we find that an in- 
crease in spatial resolution reduces this drift by an amount 
that approximately scales with the number of voxels used to 
resolve a given region. This is because the relative forces 
acting on neighboring cells decrease as the voxels are more 
closely packed together, and because the sharp density gradi- 
ents present in the ID hydrostatic profile are better resolved. 

The potential at m + 1 must be estimated because the value 
of (j)'^^^ depends on p'^^^ , which is not known until the hydro 
step has been completed. This extrapolation is one of the main 
sources of error in the code because the obtained acceleration 
is only an estimate based on the previous evolution of (/). And 
because the extrapolation is only first-order accurate in time, 
smaller ti me-ste ps do not improve the accuracy of the results 
ISpringel 2010| ). The error can be somewhat reduced by using 
a higher-order extrapolation that includes N previous time- 
steps, but our tests show that including the m-2 time- step 
only leads to a few percent reduction in error relative to the 
first-order approximation. Additionally, it imposes additional 
memory and disk overhead to save the full potential field from 
the previous N time-steps. 

For cases where the object is nearly in hydrostatic equi- 
librium the gravitational field does not rapidly change with 
time. However, a tidal disruption can lead to variations in the 
gravitational field of order unity over a fraction of the planet's 
dynamical timescale. As a result, the default method used by 
FLASH to apply the gravitational field can lead to substantial 
changes in the total energy. Motivated by this, we implement 
a novel method where the potential contribution from the fluid 
and from tidal forces are separated, negating much of the er- 
ror associated with the extrapolation of the potential. Details 
of this algorithm are provided in Appendix [A| 

3. SIMULATION RESULTS 
3.1. Single Passage Encounters 
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Figure 3. Slices through the orbital plane shortly after each periastron pas- 
sage for the simulation where rp = 2.1 rt. All plots show log p. The up- 
per, red color-coded figures show a wide view of each encounter, with white 
corresponding to p = 10~^ g cm~^ and black corresponding to p = 10~^ g 
cm~^, while the lower, blue color-coded figures show a close-up view of the 
core, with white corresponding to the maximum density pmax and black cor- 
responding to p= 10~^ g cm~^. 
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Figure 4. Same as Figure|3[ but with frames 2-10 corresponding to apastron. 
The first frame shows the initial conditions. 



For comparison purposes, we first ran a suite of simulations 
with physical initial conditions identical to that of FRW The 
planet is assumed to have a radius Rp = Rj and mass Mp = Mj, 
where Rj and Mj are the radius and mass of Jupiter. The 
planets are disrupted by a star with = lO^Mj, with the or- 
bits of incoming planets are set to have an apastron separa- 
tion Ta = 10^7?j. Our lowest-resolution models have maximum 
spatial resolutions of ^ = 0.027?j, where s is the width of the 
smallest grid cells, corresponding toN^^ 10^, slightly bet- 
ter than the peak spatial resolution of |FRW| Our results agree 
very well with FRW's results for periastron passage distances 
Tp < 2.5rt (Figure [Tj). 

Because the amount of energy stored in the oscillations is 
AE ~ E](AR/R])^ where E] is the binding energy of Jupiter 
and AR is the amplitude of the mode, the simulations con- 
verge to the analytical solution only when the spatial resolu 
tion is fine enough to resolve AR. As is evident in both FRW 
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Figure 5. Fallback accretion stream formed as the result of the 8* encounter 
between a 1 Mq star and 1 Mj planet with initial rp = 2.7rt. The large, red 
color-coded image shows a wide-view of the disruption, while the inset blue 
color-coded figure shows a close-up of the surviving planetary core. 

and our lower-resolution runs, the measured amount of en- 
ergy dissipation is larger than the analytical predictions of[PT 
when the oscillatory amplitude is smaller than the minimum 
grid scale. To test for convergent behavior, we ran higher- 
resolution simulations with double and quadrupole the res- 
olution of our fiducial test for the more grazing encounters. 
For the rp = 2.7rt and rp = 3rt runs, the improved spatial res- 
olution allows us to recover the analytical solution. It is also 
apparent that our rp = 4rt simulation is closer to convergence 
than the lower-resolution models, but the estimated resolution 
required for true convergence (/ ~ 10~^ R]_,N^ = 10^^) means 
that recovering the predicted results of|PT|for such a grazing 
passage is currently beyond our com putational ability. 

As in FRW and Khokhl ov et al .|(fl993a), we also find that 
the change in Eo^h and /spin is slightly larger than what is 
predicted by the analytical models, even in the simulations 
that have surely converged and have minimal mass loss (i.e., 
2rt < rp < 2.7rt). This is almost certainly due to the dynam- 
ical tide effects neglected by the PT model. Comparing the 
hydrodynamic results to that of a dynamical treatment of tides 
shows better a greement for this range of pericenter distances 
( Lai et al. 1994 ), but note that these dynamical models only in- 
clude the 1 = 2 /-mode of oscillation, and thus do not account 
for energy transferred to higher-order /-modes or /7-modes. 
For grazing encounters, the dynamical tide is much less im- 
portant, and thus we expect that the change in orbital energy 
and angular momentum should converge to the PT prediction 
for an inviscid planet. 

3.2. Multiple Passage Encounters 

The initial conditions of FRW are not appropriate for inves- 
tigating multiple-passage encounters as the orbital timescale 
is many thousands of times longer than the dynamical 
timescale. Thus, for our second set of simulations in which 



we explore multiple passages, we assume that the planets are 
initially on ^ = 0.9 orbits (Figure [2]). Depending on the initial 
pericenter distance, the planets are destroyed in as few as one 
and as many as ten orbits (Figures [3] and [4]), where a planet is 
considered to be destroyed once it nas lost more than 90% of 
its original mass. 

To ensure that our assumption of a less eccentric orbit is 
valid, we consider the effects of varying e when e ^ 1. For 
nearly -parabolic orbits, the shape of the orbit in the vicin- 
ity of periastron changes very little with changing e, with the 
main effect being that the average strength of the tidal force is 
slightly stronger for smaller e before and after pericenter. This 
means that the critical distance for which a planet is destroyed 
or ejected is very slightly larger than our setup for single en- 
counters described in the previous section where e o:^ 1. To 
estimate the magnitude of this effect, we calculate the ratio of 
the distances at fixed true anomaly for two orbits with eccen- 
tricities ei and e2, assuming both orbits have the same rp 



n 



ei{e2 + l)r^ 



(^i-^2)n+(^i + l)^2^p 



(1) 



As the strength of the tidal force is ex r~^ , the difference in the 
strength of the tidal force between the two orbits with respect 
to the tidal force experienced at pericenter is 



^1 



3(^1-^2) 

ei (^2 + 1) r] 



(2) 



1 +0(2) : ^1-^2^0. (3) 



This expression is maximized at ri = |rp, where the force dif- 
ference between an =0.9 and e2 = l.l encounter evaluates 
to 4%. For all other values of ri, the force differences are 
much smaller than this maximum. Thus, our hydrodynamical 
treatment of tides in a ^ = 0.9 orbit is directly applicable to all 
orbits with 0.9 <^< 1.1. 

After an orbit in which the planet sheds mass, some of the 
material that is removed from the planet's surface remains 
marginally bound to the planetary core. The majority of this 
material is then re-accreted by the planet over a few dynami- 
cal timescales. When the free-falling material encounters the 
surviving planetary remnant, its kinetic energy is converted to 
internal energy in a standing accretion shock, which results in 
the planet possessing a hot outer layer with temperature close 
to the virial temperature (Figure |5]). Additionally, the material 
striking the remnant leads to some heating of the remnant's 
outer mass shells. This effect is predominantly important in 
determining the envelope's temperature early in the accretion 
history before the pressure of the growing hot atmosphere be- 
comes comparable to the ram pressure an d is able to halt th e 
flow before it can reach the core's surface ( [Frank et al.|2002| ). 

Because the re-accreted material is marginally bound to 
planet, the orbital trajectories characterizing the accretion 
streams have apocenter distances comparable to the planet's 
Hill sphere rn = rp(Mp/3M*)^/^ and follow highly eccentric 
orbits. As the re-accreted material returns on a Keplerian t ra- 
jectory fKocha nek|1994[ |Ramirez-Ruiz & Rosswog'2009'), it 
carries a substantial amount of specific angular momentum. 
This leads to a rapid spin-up of the planetary remnant's outer 
layers. In encounters with little or no mass loss, the planets 
spin slowly post-encounter as the angular momentum carried 
by the normal modes is almost equally distributed between ro- 
tating and counter-rotating regions (Figure [6]). The process of 
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Figure 6. Pressure perturbations and spin excited by a grazing encounter 
with rp = 3rt. The upper five panels show the pressure perturbation 6P 
over several rdyn, with blue corresponding to regions of lower pressure and 
pink corresponding to regions of higher pressure relative to the base state. 
The lower five panels show the fluid's angular momentum relative to rc, 
with green representing clockwise rotation and orange representing counter- 
clockwise. Note that while the object appears to be rotating rapidly, the pres- 
ence of rotating and counter-rotating regions leads to almost an exact cancel- 
lation of the total angular momentum /. The illusion of rapid rotation is in 
fact related to the pattern speed of the / = 2, m = ±2 normal modes. Because 
the angular momentum carried by this mode is related to the tangential com- 
ponent of the expansion and contraction of the planet as it oscillates, the fluid 
vacillates back and forth at the mode frequency. As can be seen in the figure, 
the fluid with / > possess slightly more total angular momentum than the 
fluid with / < 0. This leads to a spin frequency uu = ^J/I that is actually 
very small. 



re-accretion produces a planet that has a lower average density 
and more mass at larger radii, which makes the planet easier 
to destroy on subsequent passages. 

As the re-accreted material has a temperature comparable 
to the virial temperature, radiative cooling may be able affect 
the atmosphere's structure, an effect that we do not account 
for in our simulations. However, as the planet can only ther- 
mally evolve for one orbital period, the atmosphere does not 
have much time to cool down before it has another strong tidal 
encounter with the star. While the outer layers of the planet 
may cool relatively rapidly an d lea d produce a brief transient 
visible in the UV (see Section [44| ), the denser regions of the 
planet's hot atmosphere component contains much more mass 
and is too optically thick to cool significantly before returning 
to pericenter. Assuming Thomson scattering, this corresponds 
to a mass of the hot atmosphere component Matm ~ 10~^Mj. 
This means that we expect that the thermal evolution of the 
reaccreted material is unimportant in determining the planet's 
density profile interior to the most tenuous outer layers. Thus, 
we expect that the dynamically relevant distribution of mass 
within the planet should remain relatively unchanged between 
encounters with the star, even for values of e significantly 
closer to 1 than what we use in our simulations. 

This allows us to use the change in orbital energy and an- 
gular momentum from our simulations to calculate the orbital 
evolution for orbits with semi-major axes of several au. The 
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Figure 7. Change in orbital energy £^orb attributed to each passage as a func- 
tion of Tp and orbit number A/orb- The orange regions show decreases in £^orb 
(more bound), while the cyan regions show increases in ^'orb (less bound). 
The height of each rp column shows the number of orbits a planet survives 
before being destroyed. 

change in orbital energy as a result of each encounter is shown 
in Figure [7] For grazing encounters with little mass loss, the 
change in orbital energy is negative; the planet becomes pro- 
gressively more bound to the star after each encounter. The 
magnitude of this change is related to the state of the dynam- 
ical tide on the planet at pericenter, where the interaction be- 
tween the oscillation of the fundamental modes can interact 
with the tidal field to reduce the amount of energy removed 
from the orbit, or even change its sign (see Section [33] ). As 
the mass loss per orbit exceeds ~ 10%, the trend in orbital en- 
ergy change becomes positive, and planets become progres- 
sively less bound on each subsequent encounter. This is pri- 
marily a result of the asymmetrical mass loss, although the 
interaction with the normal modes is important for encoun- 
ters where the mass loss is on the order of a few percent. 

For orbits with a ~ aice, where aice = 2.7(L*/L0)^/^ au ( |lda| 
& Lin|[2008 l) is the ice line, the total orbital energy is small 
relative to the self-binding energy of the planet, and a positive 
AE can lead to the surviving planet becoming unbound from 
the host star. The number of orbits completed by a planet be- 
fore it becomes unbound is shown in Figure [8] All Jupiter-like 
planets that come within rp < 2.7rt are expected to be ejected 
if ^ > 0.97, and destroyed otherwise. For the most grazing 
encounters and smallest eccentricities, we find that the planet 
can survive for as many as ten orbits before being destroyed 
by the host star. The general trend is that planets that come 
closer to their parent stars tend to be ejected or destroyed in 
fewer orbits, but the consequences of the encounter seem to 
depend heavily on the orientation of the planet's long axis at 
pericenter, which as we will show in the next section is quite 
difficult to predict. 

3.3. The Role of Chaos 

As the planet continues in its orbit after its first encounter 
with the parent star, the planet exhibits both rotation and os- 
cillation. Both the magnitude and sign of the orbital energy 
change is related to the phase of the planet's dynamical tide at 
periastron. Because the tidal forces strongly excite the / = 2 
modes, the coherence or decoherence of these modes at pe- 
riastron can grea tly change the effects of th at particular en- 
count er ( |Kochanek|1 992: Mardling 1995; Usami & Fujimoto 
1997| ). As the ratio of the orbital period to the break-up ro- 
tation period for a Jupiter-like planet scattered in from the 
ice line is ~ (Mj/M*)^/^(aice/^j)^^^ ~ 10^, even a very small 
change in the orbital parameters introduces a dramatic vari- 
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Figure 8. Criteria leading to planet ejection given an initial periastron pas- 
sage distance rp and eccentricity e. The colored regions correspond to the 
number of orbits before a planet is ejected by its interaction with the parent 
star. The dashed lines show the value of e given an apocenter distance equal 
to the ice line aice, assuming that L* oc M^-^. For all values of rp shown 
and for g < 0.97, planets are ejected before they are totally disrupted. Note 
that the region bounding the planets that are ejected on the first orbit has a 
monotonic dependence on rp, while all other regions exhibit more compli- 
cated structure. This is because the first passage involves a planet with no 
internal motions, and thus no relation between phase and the excitation or 
de-excitation of fundamental modes. All future passages forN> 1 involve a 
planet that is both differentially rotating and oscillating, resulting in a large 
variance in the amount of energy added or removed from the orbit. 



ability in the amount mass lost and the change in orbital en- 
ergy. 

For our simulations where e = 0.9, the planet completes 
hundreds revolutions between periastron passages, which 
means extremely fine sampling of rp would be required to 
completely describe the problem. To explore the chaotic be- 
havior we ran another multiple disruption simulation with 
Tp = 2rt, with the only difference being that the initial eccen- 
tricity is set to ^ = 0.90012, corresponding to an orbital period 
that is one free-fall timescale % = {R] /GM])^/^ /2tt longer 
than the corresponding ^ = 0.9 simulation. As shown in Figure 
|9] multiple outcomes are possible for even slight changes in 
the initial conditions, with the planet surviving for a different 
number of orbits in the two simulations. The chaotic behavior 
is also evident for our multiple disruption simulation where 
Tp = 1.8rt, which is destroyed on its fourth periastron passage, 
whereas both the rp = 1 .7rt and rp = 1 .9rt runs are destroyed on 
their third passages. This chaotic behavior arises because the 
angle of the major axis of the oscillating planetary core can 
range from being perpendicular to parallel to the true anomaly 
at periastron, which dramatically affects the ability of a planet 
to survive on any particular orbit. 

As the principle component of the tidal field are the / = 
2,m = ±2 harmonics, the only way to avoid chaotic behavior 
is to somehow remove energy from these modes before the 
planet returns to periastron. For fully-convective, Jupiter-like 
planets, this can be potentially achieved through three types of 
mechanisms: Viscosity (either microscopic or turbulent), cou- 
pling to other oscillatory modes, or the increase in entropy as- 
sociated with sound waves steepening into shocks at either the 
surface of the planet or within its interior. It should be empha- 
sized that most of the work investigating these energy-sharing 
mechanisms have concentrated on systems where the oscilla- 
tions can be treated linearly, and thus the results of these stud- 
ies can only give us a rough idea to their importance when 
applied to partially-disruptive encounters. 
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Figure 9. The evolution of two multiple encounter simulations with almost 
identical initial conditions. Both simulations use the same initial conditions 
as the rest of the multiple encounter simulations presented in this paper, with 
rp = 2rt. The only difference between the two simulations is the initial ec- 
centricity: The solid curves show the outcome for initial eccentricity e = 0.9, 
while the dotted curves show the outcome for e = 0.90012, which corresponds 
to an orbital period that is one free-fall time-scale % longer than the e = 0.9 
case. The four plots show the planetary mass Mp, orbital energy ^^orb^ spin an- 
gular momentum /p, and precession of the orbit uu as functions of t/tch- Note 
that while the outcome of the first encounter at r = is almost completely 
identical in both simulations, the behavior diverges on the second encounter, 
which leads to the planet on the ^ = 0.9 orbit being destroyed in three orbits, 
whereas the planet on the e = 0.90012 is destroyed in four. This divergence 
is a result of the phase difference between the two simulations introduced by 
the slight difference in initial orbital period. 

The two main differences between the linear models and 
the survivors of a partially-disruptive encounter are the am- 
plitude of the oscillations, which are close to unity, and the 
presence of a hot, optically thick envelope that accumulates 
after the disruption and sits on top of the oscillating core. We 
know that in order for a particular mechanism to result in sig- 
nificant damping of the / = 2 modes, the damping timescale 
has to be at least on the order of the orbital period P, if not 
shorter. As all of the mechanisms for removal of energy from 
the / = 2 modes result in a cascade to microscopic scales, sys- 
tems with effective damping will experience inflation of the 
core, inflation of the envelope, or inflation of both regions. 
This inevitably leads to reduced survivability on subsequent 
passages. 

We now briefly discuss the applicability and viability of 
each of these proposed mechanisms. For Jupiter-like plan- 
ets, the microscopic and turbulent fluid viscosities seem to 
be too small to produce any significant damping on an or- 
bital timescale ( Guil lot et al.|2 004 ). Perhaps a more promis- 
ing mechanism is the coupling of the primary 1 = 2 modes to 
higher-order "daughter" modes, which then couple to "grand- 
daughter" modes, etc., in a cascade resembling the cascade 
of energy from large scales to small scales in turbulent fluids 
(Kumar & Goodman" 1996 ). In the linear regime, the funda- 
mental mode normally couples to the low frequency g-modes, 
with the degree of coupling being related to the amplitude 
of the primary perturbation relative to the size of the object, 
AR/R. But for Jupiter-like planets, which are fully convec- 
tive and have a negligible luminosity, the polytropic index F 
is equal to the adiabatic index 7, which makes th eir interiors 
incapable of supporting g-modes (Co wling||1941| . Coupling 
can still occur through p-modes, which have higher frequen- 
cies than the fundamental mode, but the coupling is only ef- 
fective for large displacements where the behavior becomes 
non-linear a nd for which the rate of energy-sharing is highly 
uncertain ( Kum ar & Goodman] 1996 ). 

Inertial waves may be an effective means of dissipating the 
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Figure 10. Mass loss history of multiple passage encounters for different 
initial values of rp. Each curve is color-coded to correspond to a particular 
orbit number. The solid lines show the aggregate mass accreted by the star, 
while the dashed lines show the total mass lost from the planet. 



/ = 2 modes given low-frequency tidal forcing ( jivanov & Pa-| 
[paloizou |201Q| ), but the efficiency of this dissipation as the 
forcing frequency approaches the characteristic frequency is 
not well understood. And because inertial waves are most 
effective when the planet spin frequency is comparable to 
the orbital frequency, they may not be important during the 
first few passages before syncronicity is established. Addi- 
tional dissipation may occur via the interaction bet ween iner- 
tial waves and Hough waves ( Qgilvie & Lin|2Q04| ), in which 
effective coupling is achieved when the wavelength of the in- 
ertial modes is comparable to the size of the radiative zone. 
As scattered planets tend to have large eccentricities, stellar 
insolation is unlikely to be important for these planets prior to 
circularization, but as the hot envelope produced by partially 
disruptive events is radiative, effective dissipation via Hough 
waves may still be possible. 

In our simulations, we do not find that the / = 2 modes decay 
appreciably between periastron passages, despite the fact that 
the planet oscillates thousands of times in the course of each 
orbit. We do find, however, that the re-accretion of loosely- 
bound material moderates the chaotic behavior somewhat. As 
the density in this region is > 10^ times smaller than the core, 
the tidal radius for the hot envelope component is > 10 larger 
than the core's initial tidal radius, which results in the hot en- 
velope being easily removed on subsequent encounters. And 
unlike the core, the hot envelope has a large sound-crossing 
time relative to the periastron passage time, and thus no fun- 
damental modes that could affect the interaction on future pas- 
sages are excited in this region. This leads to the result that 
as the envelope becomes a larger fraction of the planet's total 
mass, the behavior becomes notably less chaotic. 

3.4. Debris Accreted by the Star 

The amount of mass removed on the first few orbits can vary 
over many orders of magnitude for a relatively small range of 
Tp. For our closest simulated encounters, nearly half the planet 
ends up being bound to the host star after the first passage, 
while our most grazing encounters show almost no mass loss 
at all (Figure 10). However, as the initial passages excite fun- 
damental modes of oscillation, and the returning debris leads 
to additional heating in the planet's outer layers, the amount of 
mass removed geometrically increases with each subsequent 
passage. This leads to the result that approximately 20-50% 
of the planet's initial mass is accreted by the host star by the 



time the planet has been completely disrupted, with slightly 
less material being accreted when the initial encounter is graz- 
ing. 

The reason closer disruptions lead to more mass accreted 
by the central star has to do with the asymmetry of the tides; 
the force resulting from the inner tide (closest to the star) at 
pericenter is 



(2rp-7?p) (7?p + rp 
(2rp+7?p) (7?p-rp 



1 + - 



37?p 



(4) 



times larger than the outer tide (furthest from the star). And 
while the effective size of the planet increases due to the ex- 
citation of oscillations and the loss of mass, subsequent pas- 
sages yield diminishing returns as material has already been 
removed on previous encounters. 

As the mass lost by the planet from a series of partially 
disruptive encounters is highly stochastic, the exact amount 
of mass accreted by the star for any given multiple-orbit en- 
counter is difficult to predict. To parameterize the average 
mass accreted by the star, we consider two limiting cases. If 
the planet is scattered into an orbit with a <C ^ice, the planet 
will be completely destroyed before it is ejected. On the other 
hand, if a planet is scattered from the ice line, £^orb is much 
smaller than the planet's binding energy and the planet will 
be ejected on the orbit for which A£'orb > ^orb- A fit of AM^ 
for these two limiting cases yields 

{1 .26 exp [-0.79/3"^] Mj : ao <C aice 
(5) 
9.62 exp [-2.59/3-^] Mj : ao - ^ice, 

for0.37</3<0.83, where the impact parameter /3 = rt/rp. 
For the case where the planet is completely destroyed, the 
amount of mass accreted is nearly constant, with a slight 
decrease with decreasing j3. The decrease is substantially 
steeper for the incomplete disruptions, as planets on grazing 
orbits are less likely to transfer much mass to their parent stars 
prior to being ejected. 

4. DISCUSSION 

4. 1 . The Jupiter Exclusion Zone 

As discussed in Section [3] there exists an exclusion zone 
with radius rr = 2.7 rt within which all Jupiter-like planets are 
either ejected or destroyed. For our simulations we used a 
polytropic model of a gas giant, with a mass and radius equal 
to present-day Jupiter. As some gas giant planets may not 
have had much time to cool before being disrupted, it is likely 
that disrupted Jupit er-like planets are larger th an our fiducial 
cold Jupiter model ( Bodenheimer et aL]|2001| ). This is con- 
firmed by observations, many of the known hot Jupiters are 
observed to be inflated, with larger radii and lower densities as 
a result of the injection of entropy via stellar inso lation ( |Fort^ 
ney et al.|2007t [Miller et al.|2009t[Lret al.|2010|), and poten 



tially also other oi s sipative mechanisms ([Qgilvie & Lin 2004 
[Laine et al^[2009t [Arras & Socrates |[2010| ). But this means 
that the tidal radii for these planets must be larger than what 
we use in our calculations, translating to planets that are more 
easily disrupted than cold gas giants. And while the core mass 
of Jupiter-like planets is uncertain, all plausible models in- 
clude cores that compose such a small fraction of the planet's 
total mass that they are unimportant in determining the struc- 
ture of the gas envelope, and thus the dynamics of the disrup- 
tion. The value of r^ calculated in this paper is thus a lower 
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Figure 11. Possible apastrons for the known hot Jupiters with Mp > O.lMj. Each arc shows the apocenter of an orbit with pericenter rp, assuming that the 
angular momentum of the initial orbit is equal to the orbital angular momentum observed today. The open circles show each planet's currently observed rp and ra, 
scaled to the tidal radius n and the ice line a\c& respectively, which are determined by the planet/host star properties. Blue-filled circles are planets that are thought 
to be aligned with their host stars, red-filled circles are thought to be misaligned (either via direct measurement or because they have observed to have statistically 
inconsistent rotation rates, see Schlaufman 2010 1, and white-filled circles show systems with unknown orientations. The thick vertical line (labeled r^- sim) shows 
the minimum possible value ot corresponding to the Jupiter exclusion zone, which is determined through our numerical simulations for Jupiter-like planets to 
be 2.1 n, while the thinner vertical line (labeled r^- obs) shows the maximum possible value for rr, which is defined by the planet that is currently observed to be 
closest to its classically defined tidal radius, WASP- 19 b. Filled black circles show the intersection between the Vr lines and each of the arcs of constant angular 
momentum shows the maximum apastron distance ra,max a planet could have been scattered from without being destroyed on its subsequent encounters with the 
host star. Note that the typical ra,max values are significantly smaller than ^ice, which indicates that those planets must have migrat ed prior to being scatt ered 
if planet-planet scattering brought them to their current positions. All data taken from the Extrasolar Planets Encyclopaedia (http : //exoplanet . eu) and 
Rene Heller's Holt-Rossiter-McLaughlin Encyclopaedia (http : / /www . aip . de/People/RHeller/ » on December 8, 2010. 



limit for Jupiter-like planets. For gas giants with masses more 
similar to Neptune, the affect of an increased core mass may 
allow these planets to survive closer to their parent stars due 
to the increase in average density. We stress that our model 
should only be applied to planets of M > O.lMj. 

Because the spatial resolution required to resolve grazing 
passages becomes computationally prohibitive beyond ~ 3rt, 
the true value of may lie beyond what we have been able 
to calculate for even cold Jupiter-like planets. Thus, the verti- 
cal cut-off line denoted by r^- sim in Figure [TT] almost certainly 
lies closer to the right of the diagram, whicn would add more 
planets to the list that could not have been directly scattered 
to their present positions from the ice line. However, the true 
cutoff value certainly lies no further than the current pericen- 
ter distance of WASP- 19 b (denoted as r^- obs). currently the 
exoplanet with the small est known value of r^/n, which ap- 



Because specific orbital angular momentum is very nearly 
conserved during even strong tidal encounters (Figure [2]), the 
currently observed semi-major axes of hot Jupiters is at most 
2 Tp^o, where rp^o is the planet's initial pericenter distance after 
being scattered into a disruptive orbit. As the exclusion radius 
Vj- is larger than the a/ 2 for many known exoplanets, there ex- 
ists a maximum initial apastron distance from which a planet 
could have been scattered without having been destroyed or 
ejected 



^obs(l-gobs)^r 



(6) 



where ^obs and ^obs are the currently observed semi-major axis 
and eccentricity. This maximum apastron value is illustrated 
for the currently known hot Jupiters in Figure [TTl Because a 



pears to be quite inflated ( [Hebb et al.J 201Q ) and probably has 
a tidal radius that is larger than a cold gas giant of the same 
mass. Hence, we can only constrain r^- to lie within the range 
^r,sim <rr< rr,obs, or 2.7 < Vr/n < 3.54. 



4.2. Implications for the Formation of Hot Jupiters 



planet that is scattered into a disruptive orbit from the ice line 
possesses ^ 1, the maximum allowed apastron is a highly 
sensitive function of r^- for r^ > O.laice. Thus, if the current 
distance of an exoplanet from its host star is less than Ir^, 
it is highly likely that its initial orbit prior to scattering must 
have been significantly closer to the star than the ice line. This 
means that no Jupiter-like planets would be observed for ini- 
tial pericenter distances smaller than this value if they were 
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Table 1 



Planet 


rp/rt 


^ a, max/ ^ice 


e* 


a 

max 


'HifCjmax^ 












(Gyr) 


HAT-P-12b 


5.2 


0.41 


1 X 10"^ - 


- 1 X 10^ 


0.4 — 4 


OGLE-TR-56 b 


4.8 


0.029 


2 X 10^ - 


-4x 10^ 


0.9 — 2 


WASP-4 b 


5.1 


0.093 


5 X 10^ - 


- 2 X 10^ 


1 — 6 


WASP- 12 b 


4.7 


0.021 


5 X 10^ - 


- 2 X 10^ 


0.5 — 1 


WASP- 19 b 


3.5 


0.010 


2 X 10^ - 


- 2 X 10^ 


0.07 — 1 



^ Assuming r^^o = 2rr, flo = '^ice- 
are calculated usin g the lower and 
|SchlaufmanH20T0t . 

scattered from the ice line to their present-day orbits. 

Five of the currently known hot Jupiter s (HAT-P-12 b, 
OGLE-TR-56 b, WASP-4 b, WASP-12 b, and WASP-19 b) 
have observed semi-major axes ^obs > 2rT- sim, corresponding 
to ra,max < ^ice (Table For four of the five planets, the 
maximum initial apastron distance is less than 0.1 aice- As 
all Jupiter-mass planets are thought to form beyond <2ice, these 
planets could not have been directly scattered to their current 
locations. This implies that a migration process that resulted 
in a dramatic decrease in a must have occurred either before 
or after these planets were scattered. If the true value of r^- 
for cold Jupiter-like planets lies closer to the maximum pos- 
sible value set by WASP-19 b (r^- obs), as many as 18 of the 
currently known hot Jupiters with M > O.lMj could not have 
survived the scattering event that deposited them at their cur- 
rent locations. 

As both planet-planet scattering and the Kozai mechanism 
do not radically alter a planet's semi-major axis, disk migra- 
tion would be required to explain how close-in planets would 
have initial apastrons that are sig nificantly smaller than aicQ 
( |Lin etal.l 19961 [Ida & Lin|20Q4||2008^ ). As many of the cur- 
rently known exoplanets have r^'max < ^ice. this would imply 
that the migration timescale in these systems must have been 
substantially shorter than the disk lifetime if the planets mi- 
grated prior to the scattering event. For the misaligned sys- 
tems, a dynamical process (either planet-planet scattering or 
the Kozai mechanism) may need to occur after the migra tion 
or while the gas disk dissipates ( [Matsumura et al.|[2Q10| ) to 
produce the observed misalignment. 

Alternatively, the planets may have migrated after the scat- 
tering event. For this to be the case, the scattering event had 
to bring the planet close enough to its parent star such that 
the tide raised on the star can alter the planet's orbit in a time 
shorter than the system age, but not too close that the planet 
is destroyed or ejected by the star. Suppose that a planet is 
scattered from beyond aicQ such that rp o > . After circular- 
ization the pericenter distance doubles to Ir^, and the planet's 
semi-major axis evolves via the interaction between the planet 



The minimum and maximum values 
upper limits for rage as compiled by 



and the tide it rai ses on the surface of the star (Eggleton et al. 
T9981 [Hut|p^80| ). At this distance, the orbital period of the 



planet is almost always shorter than the spin-pe riod of the star, 



except perhaps for stars with ages < 650 Myr ( Irwin & Bou 
lyier 2009 ). This results in a spin-up of the star, which tries to 
"catch up" to the Keplerian frequency imposed by the planet's 
orbit, and an inward migration of the planet. 

As the planet is likely to be tidally locked even prior to cir- 
cularization, the timescale for evolution of the planet's semi- 
major axis is entirely determined by the star's properties and 



the orbital frequency oo ( [Dobbs-Dixon et al.|2004| ), 

5 



a 1 
a 9 



Mp 



(7) 



where is the star's radius, 2* is the star's tidal quality fac- 
tor, and ft^ is its rotation frequency. The fastest inward mi- 
gration occurs when a star is not rotating, e.g. ^ 0. As 
planet-planet scatte ring seems to be rare after 10^ yr (iMat-J 
sumura et al.||2008| , we can sQt a/d equal to the system age 



Tage to determine (^*,max, the maximum tidal quality factor for 
which a planet can migrate from to ^obs 



: 7 X 10^ 



Mj 



8/3 



-8/3 



^0 A ^ f ^age \ 



3 days J \ Gyr J ' 



(8) 



where Pq is the initial orbital period. When setting = r^-^sim, 
all of the known hot Jupiters with ^obs < ^r-^ yields values for 
2* max that are consistent with those expected for stars (Table 

If the true radius of disruption lies closer to r^- obs. the 
measured values of 2*, max are more restrictive, meaning that 
higher dissipation rates would be required to enable planet- 
planet scattering to viably produce hot Jupiters with a < Ir^ . 
This would also imply that most observed hot Jupiters would 
only exist a short while longer before being destroyed by their 
parent star. This remaining lifetime can be calculated by solv- 
ing Equation (|7| using present-day values for a, u, and ft^ 
(denoted riife,max in Table [T]). As 2*, max is extremely sensi- 
tive to the exact value of r^-, only a slight increase in r^- is 
required to eliminate direct scattering as a possible genesis 
mechanism for a significant fraction of the hot Jupiters. If it 
can be shown that is substantially larger than what we have 
calculated in this work (i.e. — r^^ohs), 2*,max can put defini- 
tive constraints on the mechanism responsible for producing 
hot Jupiters. 

4.3. Stellar Spin-Up from Planetary Disruption 

As discussed in Section |3.4[ the star acquires a sub- 
stantial injection of angular momentum from the disrupted 
planet. Our simulations give us an empirical determination of 
AM* (/3,Mp =Mj,M* = lO^Mj), the amount of mass gained 
by a star with M* = lO^Mj from the disruption of a Jupiter- 
mass planet given the impact parameter /3 (Figure [TO]). The 
angular momentum acquired by the star is simply equal to the 
specific angular momentum at the star's radius multiplied by 
the amount of mass accreted. A/* = AM* ^/GM^R^. Note that 
this is not the same as the total angular momentum content of 
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Figure 12. Changes to the stellar spin as a result of accreting tidally disrupted planets. The curves show the cumulative probability that a star will possess a 
given angular momentum and spin inclination V^* after 1 (solid) or 5 (dashed) planetary disruptions for different initial stellar angular momentum /*,init. 
Planets are assumed to have a logarithmic distribution in mass (O.lMj <M< lOMj) and semi-major axes (jr <a< lO^ice), where Vr is the minimum distance 
for which a planet won't be tidally disrupted and a\c& is the ice line. The eccentricity e and inclination / relative to the invariable plane are assumed to follow 
Rayleigh distributions, with (jg = 0.3 and cr, = 10°. 



the disk formed from the debris as claimed by Jackson et al. 
( [2009 ), as the accretion stream does not intersect t he star's ra^ 
dius u nless the planet's original orbit has < dKochanek 

[1994]). 

If all disrupted Jupiter-like planets are approximately the 
same size and ^ ~ 1 , all disruptions can be treated using our 
fiducial model and the inclusion of the simple scaling associ- 
ated with the increase in the orbital angular momentum. This 
is because the only difference that arises from changing the 
mass of the planet and the star is the degree of asymmetry 
of the tides, which is third-order in the force expansion and 
is {R]/ Pr^f- 10~^ times smaller than the tidal force itself. 
Thus, in the event that the planet is completely destroyed by 
the star, we expect that the amount of angular momentum ac- 
quired by the star should simply scale with and Mp 



A/*(/3,Mp,M*): 



1/2 



3/2 



X A/*(/3,Mj,10^Mj). 



(9) 



However, our results show that a substantial fraction of 
planets are ejected from the system before they can be fully 
destroyed by the star. This reduces the amount of mass AM^ 
accreted by the star, and thus A/*. Therefore A/* for any 
given encounter should depend on which particular orbit a 
planet is ejected N^y This parameter depends on the orbital 
energy £'orb, and the change in orbital energy associated with 
each encounter, which as we explain in Section |3.2| should 
also should scale simply with the masses of the planet and the 
star. Thus our expression for A/* becomes 



A/*(/3,Mp,M*) = 



1/2 



3/2 



Mp 

X A/*(/3,Mj,10^Mj,A^ei) 



(10) 

(11) 



where A/ej is equal to the lowest value of N for which Equation 



( [TT) is satisfied. While this gives us the amount of angular 
momentum acquired by the star for a set of encounter param- 
eters, we must know the distribution of (3 in order to estimate 
the probability for which a star will possess a total angular 
momentum and obliquity tp^ after a number of planetary 



disruptions A/d. 

The integrated rate of scattering (i.e., all encounters with 
/3 greater than some val ue) has be en evaluated numerically 
by a number of authors. Juric & Tremaine (2008} show that 
up to 20% of planets present at the end of phase I of plan- 
etary formatio n can collide with the host star, e.g. < R^. 
[Ford & Rasio| p008 ) show that up to 16% of planets in a 3- 
body system can be thrown into a "star-grazing" orbit, which 
they define as < 10~^ainit. Nagasawa et al.^ ( 200 8 ) show 
that planet-planet scattering events tend to induce transitions 
between Kozai states (with the duration of each state being 
~ 10^ - 10^ yr) until a planet is ejected from the system or un- 
til the eccentricity of the innermost planet is damped by tides. 
Nagasawa et al . note that most of the close-in planets seem 
to be driven to their closest approaches via the Kozai mech- 
anism, which gently drives the planets into the region where 
they can be circularized, as opposed to being directly scat- 
tered into such orbits. None of the models presented above 
include the precession associated with general relativity that 
would normally quench the Kozai mechanism. 

Because the numerical scattering experiments do not in- 
clude a hydrodynamical treatment of tides, the fates of planets 
that are either scattered or driven by the Kozai to rt < rp < 
are not accurately represented. What the scattering experi- 
ments do reveal is the total rate of planet-planet interactions 
as a function of the number of gravitating bodies in the sys- 
tem, and the distribution of orbits that arises from numerous 
planet-planet interactions. And despite the simplistic treat- 
ment of tides or the neglect of tides altogether, different pre- 
scriptions of tidal dissipation do not seem to strongly affect 
the distribution of the remaining planets in the system ( ^Na-| 
gasawa et al.|2008| ), which means that the orbit distributions 



these models predict should still be appropriate to use as in- 
puts for our post-disruption stellar spin estimates. 

For systems that are not dynamically stable after the gas 
disk dissipates, or for sys tems which are driven to instability 
by an external perturber ("Zhou et al."2007^, 'Malmberg e t al. 
2010]), the models indicate that the eccentricity distribution 



of plane ts quickly evolves to a Rayleigh distribution ( |Juric & 
|Tremaine|2008) ) 



dN = 



2 \ O 2 1 



I de 



(12) 
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As the relaxation to this distribution is mainly driven by strong 
two-body planet-planet interactions, the new eccentricity e of 
a planet after a scattering event should also be drawn from the 
above distribution. Because we are interested in objects that 
may be unbound from the system after the encounter, objects 
with e initially larger than 1 should be included when calcu- 
lating the number of events at each e. Using our definition of 

(13) 



/3nAa + i' 



and by making a change of variable from e to /3, Equation ( 12 ( 
becomes 



"2(rt + ra/3)V2 



dp. 



(14) 



When r^, as is the case in planet disruption, the above 
expression simplifies to 



dN'. 



2rt 



rexp 



1 

"2^ 



dp 



(15) 



which when integrated yields a di sruption probability that 
scales inversely with {3 ( |Rees|1988| ). Because we are includ- 
ing hyperbolic encounters, can assume both positive and 
negative values, with rt/ra < < oo. To first order, the value 
of the integrand is equal for both positive and negative values 
of /3, and the total number of events where /3 > /3Ms given by 



("^ dN 



d(3 



(16) 



This implies that equal numbers of planets will be scattered 
into prograde and retrograde orbits. Under these assumptions 
and given our empirically determined lower limit for ejection 
Tp = 2.7rt, it is immediately clear that the rate of collisions 
with the central star is lower than the combined rate of ejec- 
tions and disruptions by at least a factor r^-//?* - 1 = 1.78, as- 
suming solar and Jovian values. 

Now that we have a model for the expected initial distribu- 
tion of giant planets in a dynamically relaxed system, we can 
use Equation ( 14 ) to evaluate the expected values of and ip^ 
after A/d planetary disruptions for a star of = Mq (Figure 
T2| ). Here we consider the scattered objects to be cold. Jupiter- 
like planets with Mp > O.lMj, resulting in a nearly-constant 
planetary radius Rp. This assumption is only valid if the sys- 
tem is older than ~ 10^ yr and if the planet's initial orbit is far 
enough from its parent star to have negligible insolation, but 
note that including these effects would only act to increase the 
amount of mass removed from the planet per orbit. We also 
assume that the planets are distributed uniformly in log a from 
rr to lO^ice, and in log Mp from 0.1 to 10 Mj, with Rayleigh 
distributions in e and / with ae = 03 and a/ = 10°. For sim- 
plicity, we assume that any angular momentum acquired by 
the star through a disruption is shared equally with all parts of 
the star. 

For a single planet disruption in a system where the star 
possesses initial angular momentum equal to the Sun, i/j^ ex- 
ceeds 30° in 15% of the stellar hosts, and 90° (i.e. the star ro- 
tates retrograde to the invariable plane) in 8% of systems. The 
spin rate of the star also tends to increase, with 10% of stars 
possessing 3Jq after the disruption. If disruptions are very 
common (A/^ = 5), the probability of > 10/ > 90° increases to 
47/22%, and 41% of stars have /* > 3Jq. The probability for 



enhanced values of and i/j^ are all slightly smaller if we re- 
strict disruptions to come from a > ^ice, as the amount of mass 
the star acquires from disruptions is slightly less on average. 
Equation ([5]). However, the effect is minor, with changes in 
the cumulative probabilities being on the order of a few per- 
cent. 

If a star is unable to share the angular momentum deposited 
in its outer layers in a time less than its age, the star may 
be observed to have larger values of or i/j^ than what 
would be expected given complete mixing. The timescale 
Tj^ for sharing angular momentum across the t achocline in 
the Sun is known to be only ~3 Myr (Gough & Mclntyre 
,1998^ ), with the timescale decreasing for increasing rotation 
rates in the convective region. The timescale for sharing angu- 
lar momentum across the tachocline appears to increase sig- 
nificantly as the size of the co nv ective zone shrinks f or rapidly 
rotating stars (B ames| |2003 : Barn es & Kim||2010| ), but this 
may be moderated by a fingering instability made possible 
by the larger molecular w eight of the disruption debris ( |Ga-| 
|raud|2010}|Rosenblum et al.|2010 ). A disruption in a system 
where the host star never has a thick convective zone could ef- 
fectively erase the star's original spin rate and obliquity, with 
tp^ being pulled from the inclination distribution of planetary 
orbits. 

4.4. Observational Signatures 

If we assume that the pressure profile of the hot atmosphere 
that accumulates on the planet has a pressure profile P = Pram 
at all radii (Frank et al. 2002), the initial Kelvin-Helmholtz 
cooling timescale given a total atmosphere mass Matm is ap- 
proximately 



tkh - 1.2 



Mat 



O.lMj 



Mp 
Mj 

^atm 

3^ 



105 K 



-4 



days, (17) 



where Tvir is the virial temperature, Matm is the mass of the 
atmosphere and R^tm is the radius of the atmosphere. Because 
Tvir is initially very large, the atmosphere is at first completely 
ionized. At this temperature, tkh is a few days, leading to 
a rapid thermal evolution of the planet's outer layers shortly 
after most of the sundered mass returns to the planet. During 
this phase, the planet can briefly outshine its parent star with 
Lboi ~ 10^^ erg s~\ with most of the radiation being emitted 
in the UV. As the atmosphere cools it shrinks back down onto 
the planet's surface until 7?atm ^ ^p- We can then estimate the 
temperature at which the planet radiates for the majority of its 
orbit by setting P equal to tkh, which yields a temperature of 
a few 10^ K. This indicates that the planet's outer layers will 
still be somewhat inflated before the planet returns to pericen- 
ter, and thus will be easily removed on subsequent passages. 

For planets that are ejected from their host stars, the thermal 
evolution of their outer layers continues until the temperature 
reaches a few thousand Kelvin, at which point hydrogen be- 
gins to recombine, which acts as a thermostat to maintain the 
temperature at a relatively constant value. The recombination 
timescale is 



150 



0.7 



Ma, 



O.IMt 



years, 



(18) 



where Xn is the Hydrogen fraction. As a result of this rel- 
atively long time-scale, these ejected planets could remain 
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quite bright (Lboi ~ OALq) for an extended period of time, 
even without additional tidal forcing. If we assume that one 
in ten planet-hosting stars ejects a Jupiter-like planet via a par- 
tial disruption, and adopting an average star formation rate of 
IMq per year, at least one ejected planet in the recombination 
phase should be visible in the galaxy at any one time. 

5. CONCLUSIONS 
5.1. Limitations and Future Directions 

The principle assumptions that we have made in this pa- 
per is that Jupiter-like planets are represented accurately by a 
polytropic model of its structure. One advantage of this model 
is that disruption simulations are trivially scalable to planets 
of a different size by a simple correction to (3, assuming that 
the planet's mass interior to a given radius Mp(< r) scales self- 
similarly and that the fluid 7 remains unchanged. Ann = \ 
polytrope reproduces the mass profile of coreless 1 Mj planet 
relatively well, with the difference in Mp(< r) never exceed- 
ing 10% throughout the planet's interior (N. Miller, private 
communication) . 

The inclusion of a core of a few tens of affects Mp(< 
r) out to a few times the core radius, for which Mp(< r) ~ 
0.4M(/?p). Beyond this radius, the structure of the planet is 
nearly identical to the coreless/poly tropic models. This means 
that our simulations should be an accurate representation for 
disruptions where < 70% of the planet's mass is removed for 
Jupiter-like planets. For Neptune-like planets, where the core 
mass can be larger than the gas mass, the difference in Mp(< 
r) is substantial all the way to the planet's outermost layers, 
and thus our simulation results should not be directly applied. 
As the average densities of Neptune-like planets is larger than 
Jupiter-like planets, r^ for Neptunes should assume a smaller 
value. 

Additionally, we assume that 7 = 2 throughout the simu- 
lation volume, even for regions of very low density where 
the fluid is completely ionized and should behave as an ideal 
gas (7 = 5/3) or even a radiation pressure dominated fluid 
(7 = 4/3) in the lowest-density regions. This transition to dif- 
ferent values of 7 should affect the structure of the hot enve- 
lope that forms from the re-accreted debris that surrounds a 
partially disrupted planet, which is dynamically unimportant 
but may affect the planet's observable signature. This is not 
to say that a more realistic equation of state would not affect 
the mass loss itself. As the process of ripping material from 
the planet involves rapid fluid decompression, a decrease in 7 
may result in slightly altered disruption dynamics. 

Ideally, one would like to extend the models we have pre- 
sented here to include a more physical equation of state that 
can treat all components of the pre- and post-disrupted planet 
realistically. As the resolution required to determine r^ for 
multiple-orbit encounters beyond what we have presented 
here is prohibitive, it seems that the exploring the affects of 
using a more-complete equation of state with a realistic initial 
planet model is the next natural step for future studies. In the 
case of planets with a substantial core, these modifications are 
necessary to determine r^ with any confidence. 

5.2. The Fates of Scattered Jupiters 

The fate of a Jupiter-like planet after a strong scattering 
event is a function of the strength of the tidal forces it ex- 
periences at periastron. In this paper, we have determined 
the disruption radius r^ for Jupiter-like planets which sets the 
boundary between long-term survival and rapid tidal disrup- 



tion. Below, we summarize the various post-scattering out- 
comes in order of decreasing distance, using r^- and the tidal 
radius r^ as points of reference. 

Stalled (rp > 6rt) : The planet is deposited into an orbit 
where the rate of tidal dissipation is too small to result in a 
change in semi-major axis over the lifetime of the system. 
This planet may be in a Kozai state driven by a third body 
in the system, or could experience another strong scattering 
event, which may lead to an increase of eccentricity and sub- 
sequent circularization. 

Circularization/Migration [r^ < r^ < Sr^^e < 0.9): In this 
region, the planet is close enough to its parent star that tidal 
dissipation is effective, and the planet can circularize in 10^ yr 
or less for moderate values of e. For near-radial orbits, circu- 
larization may still be longer than the stellar age, but again the 
Kozai mechanism or scattering could lead to a more rapid or- 
bital evolution. All currently observed hot Jupiters are either 
stalled, in the process of circularizing, or already have circu- 
lar orbits. If the planet is close enough to its parent star and 
G* < 10^, the planet will raise a tide on the star and migrate 
inwards due to the transfer of angular momentum. 

Ejection (^p < r^-,^ > 0.97): A planet that passes within 
the exclusion zone will be ejected from the system if its ini- 
tial orbit is radial enough such that its orbital energy £^orb is 
significantly smaller than the self-binding energy of the 
planet. Slightly less than half of the planet's initial mass re- 
mains bound to the central star, carrying with it a large reser- 
voir of angular momentum that can significantly alter the host 
star's spin rate and axis of rotation. The ejected planet will 
remain as bright or brighter than its host star for a few years, 
eventually plateauing via hydrogen recombination as an ob- 
ject with a fraction of solar luminosity for a century. All 
Jupiter-like planets that scatter in from beyond ^ice such that 
Tp < rr will be ejected from the system. 

Complete disruption {r^ Kr^-^e < 0.97): For planets that 
are deep within their parent star's potential well, the planet 
cannot soak the change in energy required to significantly al- 
ter the orbit, which eventually leads to its complete disruption. 
Approximately half of the planet's mass accretes onto the stel- 
lar host, carrying the same specific angular momentum as in 
the ejection case, leading to even more pronounced effects 
on the stellar spin. Only planets that have migrated close to 
their stars prior to being scattered are destroyed before they 
are ejected. 

Collision with the central star (rp < 7?*) : The planet strikes 
the surface of the star directly. Anywhere from half to all 
of the planet's mass is absorbed by the star, with the angu- 
lar momentum being carried by this material potentially be- 
ing smaller than that carried by the debris from a disruption, 
depending on how direct the impact is. These events are ap- 
proximately twice as uncommon as ejections/disruptions. 

We have benefited from many useful discussions with E. 
Ford, K. Schlaufman, E. Quataert, F. Rasio, G. Laughlin, N. 
Miller, D. Fabrycky, B. Hansen, M. Rees, and A. Socrates. 
The software used in this work was in part developed by 
the DOE-supported ASCI/Alliance Center for Astrophysical 
Thermonuclear Flashes at the University of Chicago. Com- 
putations were performed on the Pleiades and Laozi UCSC 
computer clusters. We acknowledge support from the David 
and Lucille Packard Foundation (JG and ER-R), NSF grants 
PHY-0503584 (ER-R) and AST-0908807 (DL), NASA grants 
NNX07A-L13G (DL), NNX07AI88G (DL), NNX08AL41G 
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Figure 13. Cartoon showing the sequence of states used to compute the evolution over a single a time-step using our modified gravity solver. The initial state m 
is showing in the left panel, the middle panel shows the intermediate state (m-i- 1), and the last panel shows the final state m+l. r* is the position of the point mass 
representing the star, Yc is the planetary core's true center of mass, rx = (Mpr^ -i-M*r*)/(Mj -I-M*) is the center of mass of the complete system. The variable 
vector field U(r) = [p,pv,e,X] represents the value of all the conserved quantities in a given state. Differences between the states are exaggerated for illustrative 
purposes. 



(DL), and NNX08AM84G (DL); and the NASA Earth and Space Science Fellowship (JG). 

APPENDIX 
MODIFIED PPM GRAVITY ALGORITHM 

Because the binding energy of a planet on a disruptive orbit is comparable to the planet self-binding energy, the conservative 
properties of a code used to investigate planetary disruption are important. As the simulation of a partially-disruptive encounter 
involves the simultaneous resolution of both a compact core and two debris tails which are hundreds of time larger than the core, 
we found that the standard methods used to calculate the gravitational potential in a tidal disruption are too computationally 
expensive given the required accuracies. 

Our approach was to improve upon the gravity solver found within the FLASH hydrodynamics code ( [Fryxell et al.|200Q| such 
that it is better suited to investigating the problem of tidal disruption, a pictorial representation of the algorithm described below 
is shown in Figure [13] Two centers of mass are calculated at each time-step, which are defined as 

r.= ^ (Al) 

(A2) 



Z^/,p(0>/Pmax^^" 

where is the "true" center of mass and is the "core" center of mass, which only includes matter above a density cut-off 
P > /Pmax, where / is set to 0.1. The planet's virtual particle is fixed to spatially coincide with the y^ vector at all times. The 
planet's self-gravity is calculated using a multipole expansion of the planet's mass about y^ instead of Ft, which allows us to better 
approximate the planet's potential using less terms in the multipole expansion. 

Our algorithm is particularly well- suited for investigating tidal disruptions. If the expansion were performed about Ft after 
large tidal tails have formed within the simulation, y^ would be mostly determined by the position of material that is far away 
from the core, and the region with the best resolution of the potential may lie in empty space. We can estimate the number of 
terms required to represent the core's potential if the expansion is carried out about Ft instead of Yc. Assume the core has a radius 
Tc, and that the core lies a distance d = |r-rc| from the true center of mass. The angular scale of a lobe corresponding to a 
spherical harmonic of degree / is simply tt//, which means that even a first-order approximation of the core's potential requires 
an expansion with /max > 7rJ/rc. Assuming ~ 10% of a planet's mass is lost during an encounter and this material lies an average 
distance ~ 10^ Rj from the core at apocenter, d is on order 1007?j, meaning that the multipole expansion must be carried out to 
/ > 300. This is highly impractical, and thus it is much more efficient to carry out the multipole expansion about the planet's 
core, whose position is associated with the densest material in the simulation and is where the potential gradients are largest. 

The potential (j) calculated from the multipole expansion about Yc is used both to apply forces to the fluid in the simulation 
domain and to the virtual star and planet particles. A consequence of not expanding about the true center of mass is that there 
exists a non-zero force that is applied to the core. These forces are associated with the odd-/ multipole terms that usually cancel 
when the expansion is carried out about the true center of mass. Our multipole expansion does not discard these odd-/ terms, 
which allows us to confidently represent the fluid's potential using an expansion of relatively low order. 

In the FLASH code's split PPM formalism, the equations used for coupling hydrodynamics and the gravitational field for a cell 
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Figure 14. Accumulation of numerical error in the total energy E and angular momentum / for the simulations ran to compare with the results of|FRW| The 
solid lines show simulations for which s = 0.02R], the dashed lines show simulations where s = O.OIR], and the dash-dotted lines show simulations where 
^ = 5 X 10~^Rj. Note that the rate of error accumulation is greatest shortly after periastron when the planet's self-potential is rapidly varying as a function of time. 
As the planet returns to a state of quasi-equilibrium, the rate of error accumulation asymptotes to the original rate. 



/ along each of the three cartesian directions are ( [Bryan et al.|1995| ) 



At"" 



The acceleration g^^^ is calculated by extrapolating ~^ and to obtain an estimate for (j)^ 



This is in turn used to calculate 



2Axi 



7m+l _ im 



Im+1 



1 + 



im+l 



12 



im-l 



Ar 



7m+l 



-24>1 



m+l 



im+l 
^i-1 



^Pi 



where 6pi is defined as 



Spi = min -pt^,2\pT- Pti \,2\pT-pZi\)^ sign (pSi " pT-i] 



(A3) 
(A4) 

(A5) 

(A6) 
(A7) 



to enforce monotonicity in p. 

Because the potential of the star in our simulations is approximated by an analytical expression (for our simulations, a monopole 
potential), we can implement the following modification such that the component of the acceleration attributed to the star can be 
calculated to much higher accuracy. We first make the assumption that the matter distribution remains fixed over the course of 
the time-step, and then integrate the virtual particle positions forward in time from m to m+ 1. This allows us to calculate 



GM,p^, 

rm+l I 



(A8) 



an estimate of the star's contribution to the potential based on the star's approximate final position. This estimate should be much 
closer to the true value than simple extrapolation as the poten tial a t a given location has a non-trivial time-dependence. Splitting 
into two components (j)^^ (star) and (planet). Equation (|A5|) becomes 



Im+l 



1 + 



Ar 



At" 



Additionally, a correction must be made to Equation (|A6| 



m+l _ m+l . 

3i,cor 5 1 



m+l 



7 m+l 



(A9) 



(AlO) 



where g^'^^ is the acceleration experienced by the core due to the presence of odd-/ terms in the multipole expansion. Using g 
instead of a hydrodynamical step is performed according to Equations (A3 ) and ( A4), which yields pj""^^ and thus the true 



m+l 
/,cor 



contribution at m + 1 of the planet to the potential, (p^y . The positions of the virtual particles are then re-integrated from m to 
m+l, but this time using a linear interpolation of the time-evolving potential over the time-step. 

As mentioned previously, a complication introduced by this method is that the net force applied to the point about which the 
multipole expansion is carried out is non-zero. This is because we do not cancel out the acceleration applied to the extended 
object's true center of mass rt, but rather the center of mass defined by the densest material, Fc. While our correction to the 
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gravitational acceleration (Equation (AlO)) mitigates the problem somewhat, the mass density within the simulation domain 
varies as a function of time, and thus the total mass that satisfies the density criteria to be included when calculating also 
changes as a function of time. This leads to small, but measurable, changes in the reference point over the course of a simulation 
which are not entirely corrected by simply subtracting off the force applied at Fc. To ensure that angular momentum and energy 
are conserved, a correction is made to virtual particle positions over the course each time- step when re-integrating from m to 
m+ 1, 

t — t^ 

x^,,At) = xAt)-j;^^^K'-r7). (All) 

In practice, and r^"^^ are very nearly equal, with a substantial difference only arising when the extended object is close to being 
destroyed. Typical corrections are significantly smaller than the size of individual grid cells in the most highly-resolved regions. 



The accumulation of error for our single passage encounters is shown in Figure 14 In the deepest encounters where the self- 
potential changes rapidly, we find that the rate of relative error accumulation in totaTenergy and angular momentum is no larger 
than 10""^ per dynamical time for our simulations with the lowest maximum-resolution, with the planet's initial diameter being 
resolved by 50 grid cells at ^ = 0. Most of this error arises from the ejected tidal tails, which are resolved at lower resolution 
out of necessity because of the large volume they occupy. For encounters in which all the matter within the simulation is always 
resolved at highest resolution (i.e. those that do not have significant mass loss), the relative error accumulation is only 10~^ per 
dynamical time. 

We also tested the error accumulation for two sets of simulations with higher maximum resolutions, with the planet's initial 
diameter being resolved by 100 and 200 grid cells. These simulations show an accumulation of fractional error of ~ 10~^ and 
10~^ per dynamical time, respectively. 

REFERENCES 



Antonini, P., Lombardi, J., & Merritt, D. 2010, eprint arXiv: 1008.5369 

Arras, P. & Socrates, A. 2010, ApJ, 714, 1 

Barnes, S. A. 2003, ApJ, 586, 464 

Barnes, S. A. & Kim, Y.-C. 2010, ApJ, 721, 675 

Bicknell, G. V. & Gingold, R. A. 1983, ApJ, 273, 749 

Bodenheimer, P, Lin, D. N. C, & Mardling, R. A. 2001, ApJ, 548, 466 

Bryan, G. L., Norman, M. L., Stone, J. M., Cen, R., & Ostriker, J. P. 1995, Computer Physics Communications, 89, 149 
Carter, B. & Luminet, J. 1983, A&A, 121, 97 
Carter, B. & Luminet, J. P 1985, MNRAS, 212, 23 
Cowling, T. G. 1941, MNRAS, 101, 367 

Diener, P, Frolov, V. P, Khokhlov, A. M., Novikov, 1. D., & Pethick, C. J. 1997, ApJ, 479, 164 

Dobbs-Dixon, L, Lin, D. N. C, & Mardling, R. A. 2004, ApJ, 610, 464 

Eggleton, P P, Kiseleva, L. G., & Hut, P 1998, ApJ, 499, 853 

Evans, C. R. & Kochanek, C. S. 1989, ApJ, 346, L13 

Faber, J. A., Rasio, F A., & Willems, B. 2005, Icarus, 175, 248 

Fabrycky, D. & Tremaine, S. 2007, ApJ, 669, 1298 

Ford, E. B. & Rasio, F A. 2008, ApJ, 686, 621 

Fortney, J. J., Marley, M. S., & Barnes, J. W. 2007, ApJ, 659, 1661 

Foucart, F & Lai, D. 2010, arXiv 

Frank, J., King, A., & Raine, D. J. 2002, Accretion Power in Astrophysics 

Frolov, V. P, Khokhlov, A. M., Novikov, 1. D., & Pethick, C. J. 1994, ApJ, 432, 680 

Fryxell, B., Olson, K., Ricker, P, Timmes, F X., Zingale, M., Lamb, D. Q., MacNeice, P, Rosner, R., Truran, J. W., & Tufo, H. 2000, ApJ, 131, 273 
Garaud, P 2010, arXiv 

Gough, D. O. & Mclntyre, M. E. 1998, Nature, 394, 755 

Guillochon, J., Ramirez-Ruiz, E., Rosswog, S., & Kasen, D. 2009, ApJ, 705, 844 

Guillot, T., Stevenson, D. J., Hubbard, W. B., & Saumon, D. 2004, In: Jupiter. The planet, 35 

Hebb, L., CoUier-Cameron, A., Triaud, A., Lister, T., Smalley, B., Maxted, P., Hellier, C, Anderson, D., PoUacco, D., Gillon, M., Queloz, D., West, R., Bentley, 

S., Enoch, B., Haswell, C, Home, K., Mayor, M., Pepe, F, Segransan, D., Skillen, I., Udry, S., & Wheatley, P 2010, ApJ, 708, 224 
Hubbard, W. B. 1984, New York 
Hut, P 1980, A&A, 92, 167 
Ida, S. & Lin, D. N. C. 2004, ApJ, 604, 388 
— . 2008, ApJ, 685, 584 

Irwin, J. & Bouvier, J. 2009, lAU Symposium, 258, 363 
Ivanov, P B. & Novikov, I. D. 2001, ApJ, 549, 467 
Ivanov, P B. & Papaloizou, J. C. B. 2010, MNRAS, 407, 1609 
Jackson, B., Barnes, R., & Greenberg, R. 2009, ApJ, 698, 1357 
Juric, M. & Tremaine, S. 2008, ApJ, 686, 603 

Khokhlov, A., Novikov, I. D., & Pethick, C. J. 1993a, Astrophysical Journal v.418, 418, 181 
— . 1993b, Astrophysical Journal v.418, 418, 163 

Kobayashi, S., Laguna, P, Phinney, E. S., & Meszaros, P 2004, ApJ, 615, 855 

Kochanek, C. S. 1992, ApJ, 385, 604 

— . 1994, ApJ, 422, 508 

Kozai, Y. 1962, AJ, 67, 591 

Kumar, P & Goodman, J. 1996, ApJ, 466, 946 

Lai, D., Rasio, F A., & Shapiro, S. L. 1994, ApJ, 437, 742 

Laine, R., Lin, D., & Dong, S. 2009, ApJ, 685, 521 

Lee, W. H., Ramirez-Ruiz, E., & van de Ven, G. 2010, ApJ, 720, 953 

Li, S.-L., Miller, N., Lin, D. N. C, & Fortney, J. J. 2010, Nature, 463, 1054 



Ejection and Disruption of Giant Planets 



17 



Lin, D. N. C, Bodenheimer, P., & Richardson, D. C. 1996, Nature, 380, 606 
Lodato, G., King, A. R., & Pringle, J. E. 2009, MNRAS, 392, 332 
Loren-Aguilar, R, Isern, J., & Garcia-Berro, E. 2009, A&A, 500, 1193 
Malmberg, D., Davies, M. B., & Heggie, D. C. 2010, arXiv, astro-ph.EP 
Mardling, R. A. 1995, ApJ, 450, 722 

Matsumura, S., Takeda, G., & Rasio, F. A. 2008, ApJ, 686, L29 

Matsumura, S., Thommes, E. W., Chatterjee, S., & Rasio, R A. 2010, ApJ, 714, 194 

Miller, N., Fortney, J. J., & Jackson, B. 2009, ApJ, 702, 1413 

Nagasawa, M., Ida, S., & Bessho, T. 2008, ApJ, 678, 498 

Nolthenius, R. A. & Katz, J. I. 1982, Astrophysical Journal, 263, 377, a&AA ID. AAA032.066.185 
Ogilvie, G. L & Lin, D. N. C. 2004, ApJ, 610, 477 

Press, W. H., Flannery, B. P., & Teukolsky, S. A. 1986, Cambridge: University Press 

Press, W. H. & Teukolsky, S. A. 1977, ApJ, 213, 183 

Ramirez-Ruiz, E. & Rosswog, S. 2009, ApJ, 697, L77 

Rathore, Y. 2005, Proquest Dissertations And Theses 2005. Section 0037, 23 

Rees, M. J. 1988, Nature, 333, 523 

Robertson, B. E., Kravtsov, A. V., Gnedin, N. Y., Abel, T., & Rudd, D. H. 2009, MNRAS, 1774 

Rosenblum, E., Garaud, R, Traxler, A., & Stellmach, S. 2010, arXiv 

Rosswog, S., Ramirez-Ruiz, E., & Hix, W. R. 2008, ApJ, 679, 1385 

— . 2009, ApJ, 695, 404 

Schlaufman, K. C. 2010, ApJ, 719, 602 

Shen, Y & Turner, E. L. 2008, ApJ, 685, 553 

Springel, V. 2010, MNRAS, 401, 791 

Takeda, G. & Rasio, F. A. 2005, ApJ, 627, 1001 

Tasker, E. J., Brunino, R., Mitchell, N. L., Michielsen, D., Hopton, S., Pearce, R R., Bryan, G. L., & Theuns, T. 2008, MNRAS, 390, 1267 

Triaud, A. H. M. J., Cameron, A. C, Queloz, D., Anderson, D. R., Gillon, M., Hebb, L., Hellier, C, Loeillet, B., Maxted, R F., Mayor, M., Pepe, F., PoUacco, 

D., Segransan, D., Smalley, B., Udry, S., West, R. G., & Wheatley, R J. 2010, arXiv, astro-ph.EP 
Usami, M. & Fujimoto, M. 1997, ApJ, 487, 489 

Watson, C, Littlefair, S., Diamond, C, Cameron, A. C, Fitzsimmons, A., Simpson, E., Moulds, V., & PoUacco, D. 2010, arXiv 
Zakamska, N. L., Pan, M., & Ford, E. B. 2010, MNRAS, 1566 
Zhou, J.-L., Lin, D. N. C, & Sun, Y.-S. 2007, ApJ, 666, 423 



